Model 1

To reduce the number of parameters in the statistical model, we categorised different fertilisers into categories consisting of;

urea: Fertilisers where the main component was urea

UAN: UAN

ammonium+1: Fertilisers including a low percentage of ammonium

ammonium+2: Fertilisers with higher percentage of ammonium

The table below depicts the categorisation of fertilisers, and the number of observations in each sub-category.

In Model 1, we also examined the effects of application method (Application) on NH3 emissions. The tables below show the number of observations within each category examined.

Table below shows counts of all categories of variables:

##             
##              Broadcast Incorporated Liquid
##   urea+            128            5     16
##   UAN                0            0     31
##   ammonium+1        20            0      0
##   ammonium+2         0            0      5

We analysed the effects of Fertiliser types (urea, UAN, ammonium+1, and ammonium+2) and Application method (Broadcast, Incorporated, and Injected), using a general linear mixed effects model (lmer, R-package lme4, Bates et al. 2015). Because the observations originated from different studies, each study most often contributing more than one observation, the reference for study was used as a random effect, to account for non-independence due to repeated measurements.

We used NH3-loss (%), calculated as relative difference between ammonia emissions from a test site and a control site proportional to the control, as a response variable. There were 0 negative observations out of the total 205.

The frequency distribution of the untransformed NH3-loss measurements is shown below.

To account for non-normality, the response variable was transformed with a square root transformation. We applied a constant of 4 to allow the analysis to include negative values.

To account for heteroscedasticity, we used weighted ordinary least squares (WOLS, Olvera Astivia & Bruno 2019). We built a vector giving larger weight to observations with smaller variation. To this end, we extracted absolute residuals and fitted values from an unweighted version of Model 1 (see below for code). We then applied a regression model between the absolute residuals (response variable and fitted values (predictor variable), to estimate the relationship between the amount of variation and the data and the model estimates. To attain a vector with higher values attached to lower variation, we used the inverse of the relationship estimated in the regression model. The weights correspond to increased weight with decrease in variation, thus giving more precedence to more accurate data.

The R-code for the acquiring the weights, and the final model is shown below.

The model assumptions are checked in the figures below:

NH3-loss was different for different fertilisers (F(3, 190) =2.76, p = 0.044), as were the emissions from different application methods (F(2, 189) =3.16, p = 0.045). Overall, the measurement method did explain variation in NH3 emissions (F(NA, NA) =NA, p NA).

The pairwise-comparisons of factor-level effects can be seen below. The comparisons are made to urea (in Fertiliser) and Broadcast (in Application). The estimates are given in the sqrt(x+4) -scale.

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sqrt(NH3loss + 4) ~ Fertiliser + Application + (1 | Ref)
##    Data: NH3predict
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 664.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2268 -0.6542 -0.1125  0.5405  2.6311 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Ref      (Intercept) 0.5493   0.7412  
##  Residual             1.2806   1.1316  
## Number of obs: 205, groups:  Ref, 38
## 
## Fixed effects:
##                         Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)               4.4195     0.1731  39.5866  25.531  < 2e-16 ***
## FertiliserUAN            -0.1118     0.3826 196.7787  -0.292  0.77045    
## Fertiliserammonium+1     -0.8621     0.3087 193.7718  -2.793  0.00575 ** 
## Fertiliserammonium+2     -0.4891     0.6131 184.4702  -0.798  0.42608    
## ApplicationIncorporated  -1.1435     0.5912 197.4932  -1.934  0.05451 .  
## ApplicationLiquid        -0.6213     0.3716 181.0293  -1.672  0.09625 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) FrtUAN Frtl+1 Frtl+2 ApplcI
## FertilsrUAN  0.039                            
## Frtlsrmmn+1 -0.190 -0.005                     
## Frtlsrmmn+2  0.063  0.436  0.038              
## ApplctnIncr -0.141 -0.005  0.027 -0.009       
## ApplictnLqd -0.238 -0.730  0.204 -0.473  0.034

The effect sizes (marginal means) for each category of each factor are presented below, parametrised to urea (in Fertiliser) and Broadcast (in Application). The effect sizes are given at the back-transformed scale (NH3-loss %).

The same effect sizes (and their 95% confidence interval) are shown below:

Model 1 explains 35.57 % of the variation in the data (conditional R2), with the fixed effects explaining 7.93 % (marginal R2). The R2-values for linear mixed effects model were calculated according to Nagakawa & Schielzeth (2013).

##
# Model explanatory power
#

#
# R-squared
r.squaredGLMM(m1_lmer_p1)
##             R2m      R2c
## [1,] 0.07928294 0.355678
#
# RMSE
sqrt(sum(resid(m1_lmer_p1)^2)/insight::get_df(m1_lmer_p1))
## [1] 1.079217
#
#
# 

# Save fitted values for future models
NH3predict$M1_fit <- fitted(m1_lmer_p1)
NH3predict$M1_res <- resid(m1_lmer_p1) # For R2-calculations

Model 2 - Application rate

In model 2, as in all models following, we applied a linear model using fitted values from Model 1 as an offset variable, to account for the variation explained by the factors already estimated in Model 1. As in Model 1, we transformed the response variable (NH3-loss %) with a square root-transformation with a constant of 4. The model was also weighted the analysis with the same weights as Model 1 to account for heteroscedasticity. This same structure was used in all following models.

The sole explanatory variable in model 2 was application rate (kg/ha).

The model as it was applied can be seen below:

# Model 2 
m2_lm <- lm(sqrt(NH3loss+4) ~ Application.rate + offset(M1_fit),
           # weights=Weight,
            data=NH3predict_m2)

The model assumptions can be seen below:

Application rate explained a significant proportion of the variation in NH3-emissions (F(1, 201) =1.32, p = 0.252). There was a statistically significant positive trend between Application rate and NH3-emissions (\(\beta\)+/-SE = 0.00174 +/- 0.00152).

There are 203 observations in the analysis.

Full model output can be seen below:

## 
## Call:
## lm(formula = sqrt(NH3loss + 4) ~ Application.rate + offset(M1_fit), 
##     data = NH3predict_m2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.5285 -0.7255 -0.1033  0.5329  2.9898 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)
## (Intercept)      -0.186420   0.177382  -1.051    0.295
## Application.rate  0.001741   0.001516   1.149    0.252
## 
## Residual standard error: 1.064 on 201 degrees of freedom
## Multiple R-squared:  0.006521,   Adjusted R-squared:  0.001578 
## F-statistic: 1.319 on 1 and 201 DF,  p-value: 0.2521

Application rate explains -0.41 % of residual variation in NH3-emissions (adjusted R-squared). (Note: Adjusted R-squared is manually adjusted to incorporate the offset-variable, and thus not consistent with the adjusted R-squared from the above model output.)

Figure below is a linear approximation of the effect of Application rate on NH3-emissions (linear slope back-transformed from the square-root scale).

## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

Model 3 - Soil pH

Soil pH was binarised into categories, where soils with pH>7.0 were categorised as chalky, and ph<=7.0 as normal. The number of observations in each group and their frequency distribution can be seen below.

Model 3 followed the structure of Model 2. We added Soil pH, Fertiliser type, and the interaction between the two as fixed factors. In this model, Fertiliser types were used to identify differences between pH-categories, not to analyse the effect of Fertiliser on its own - the offset of Model 1 removes the variation explained by Fertilisers otherwise.

# Model
m3_lm <- lm(sqrt(NH3loss+4) ~ Fertiliser * pH_cat + offset(M1_fit), 
            weights=Weight,
            data=NH3predict_m3)

The heteroscedasticity and normality of errors were as below:

The effect of soil pH on NH3-emissions was dependent on the Fertiliser types (F(3, 178) =4.37, p = 0.005). There was no independent effect of soil pH on NH3 emissions (F(1, 178) =0.01, p = 0.923).

There are 186 observations in the analysis.

Full anova-table can be seen below:

Full model output can be seen below:

## 
## Call:
## lm(formula = sqrt(NH3loss + 4) ~ Fertiliser * pH_cat + offset(M1_fit), 
##     data = NH3predict_m3, weights = Weight)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3362 -0.7841 -0.1115  0.8140  4.4495 
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                        0.04187    0.12024   0.348 0.728074    
## FertiliserUAN                     -0.36135    0.34153  -1.058 0.291484    
## Fertiliserammonium+1              -0.84222    0.35546  -2.369 0.018888 *  
## Fertiliserammonium+2              -0.17150    0.69973  -0.245 0.806668    
## pH_catchalky                      -0.33912    0.18344  -1.849 0.066160 .  
## FertiliserUAN:pH_catchalky         0.83519    0.43012   1.942 0.053747 .  
## Fertiliserammonium+1:pH_catchalky  1.45799    0.43541   3.349 0.000991 ***
## Fertiliserammonium+2:pH_catchalky  0.50115    0.79221   0.633 0.527805    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.297 on 178 degrees of freedom
## Multiple R-squared:  0.07138,    Adjusted R-squared:  0.03486 
## F-statistic: 1.954 on 7 and 178 DF,  p-value: 0.06372

The model explained 3.34 % of the residual variation in NH3-loss.

Combining Model 1 and Model 3

For acquiring estimates for Soil pH in different fertiliser types at the same time, and same scale as the base-model, we added Soil pH and the interaction between pH and FertiliserType into the initial model. Please note this model has less data than Model 1, and as it includes new parameters, the estimates change slightly.

The effect sizes (marginal means) for the interaction can be seen below. The effect sizes are given at the back-transformed scale (NH3-loss %), parametrised to Fertiliser=urea+, Application=Broadcast, pH=normal.

The same effect sizes (and their 95% confidence interval) are shown below:

Model 4 - Soil clay content

In Model 4, we explored the effect of Soil clay content (%) on NH3-loss (%).

We used Soil clay content (%) as a continuous covariate and the sole explanatory variable in the model.

# Model
m4_lm <- lm(sqrt(NH3loss +4) ~ Clay + offset(M1_fit), 
            weights=Weight, 
            data=NH3predict_m4)

The model assumptions (heteroscedasticity and normality of errors):

Soil clay content (%) affected NH3-loss (%) negatively (\(\beta\)+/-SE = 0.00356 +/- 0.00388, F(1, 174) =0.84, p = 0.361).

There are 176 observations in the analysis.

Full model output can be seen below:

## 
## Call:
## lm(formula = sqrt(NH3loss + 4) ~ Clay + offset(M1_fit), data = NH3predict_m4, 
##     weights = Weight)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.5492 -0.7728 -0.0676  0.6767  4.9120 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.135601   0.146621  -0.925    0.356
## Clay         0.003555   0.003879   0.916    0.361
## 
## Residual standard error: 1.155 on 174 degrees of freedom
## Multiple R-squared:  0.004804,   Adjusted R-squared:  -0.0009156 
## F-statistic: 0.8399 on 1 and 174 DF,  p-value: 0.3607

The model explained 2.2 % of the residual variation in NH3-loss.

Figure below is a linear approximation of the effect of Soil Clay content (%) on NH3-emissions (linear slope back-transformed from the square-root scale).

Model 5 - Soil organic carbon

Soil organic carbon (g/kg) was modeled as a continuous variable. Most soils in the data were mineral (with a cutoff of 120 g/kg SOC), with only 3 organic soils in the dataset.

We used Soil organic carbon (g/kg) as a sole continuous covariate in the model.

# Model
m5_lm <- lm(sqrt(NH3loss+4) ~ SOC + offset(M1_fit), 
            weights=Weight,
            data=NH3predict_m5)

The model assumptions (heteroscedasticity and normality of errors):

There was no effect of Soil organic carbon (g/kg) on NH3 loss (F(1, 82) =1.03, p = 0.314).

There are 84 observations in the analysis.

Full model output can be seen below:

## 
## Call:
## lm(formula = sqrt(NH3loss + 4) ~ SOC + offset(M1_fit), data = NH3predict_m5, 
##     weights = Weight)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3239 -1.0363  0.1124  1.0401  3.8339 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.31172    0.21331  -1.461    0.148
## SOC          0.01122    0.01107   1.014    0.314
## 
## Residual standard error: 1.34 on 82 degrees of freedom
## Multiple R-squared:  0.01237,    Adjusted R-squared:  0.0003291 
## F-statistic: 1.027 on 1 and 82 DF,  p-value: 0.3138

Adjusted R2-value of the model was 0, and multiple R2 was 0.012373. The negative adjusted R2-value indicates that the amount of variation accounted for by Soil organic carbon is negligible. In practice, soil organic carbon explains no variation in the data.

Figure below is a linear approximation of the effect of Soil organic carbon (g/kg) on NH3-emissions (linear slope back-transformed from the square-root scale).

Model 6 - Cover

Crop cover was split into three categories: no cover (“none”), arable crops (including alfalfa, cereal, maize, rape and winter cereal), and grass. The number of observations in each group, and their NH3-loss (%) distributions can be seen below:

We added Cover as a fixed factor and a sole explanatory variable in the model.

# Model
m6_lm <- lm(sqrt(NH3loss+4) ~ Cover + offset(M1_fit), 
            weights=Weight,
            data=NH3predict_m6)

The model assumptions (heteroscedasticity and normality of errors):

The cover categories (none, arable, grass) did not significantly differ from one another in NH3-emissions (F(3, 201) =3.58, p = 0.015).

Full model output can be seen below:

## 
## Call:
## lm(formula = sqrt(NH3loss + 4) ~ Cover + offset(M1_fit), data = NH3predict_m6, 
##     weights = Weight)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3593 -0.8537 -0.1245  0.8459  4.7474 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   0.1405     0.1139   1.234  0.21874   
## CoverArable  -0.4019     0.1468  -2.737  0.00676 **
## CoverGrass    0.1045     0.2155   0.485  0.62837   
## CoverOther   -0.5963     1.1028  -0.541  0.58932   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.248 on 201 degrees of freedom
## Multiple R-squared:  0.05066,    Adjusted R-squared:  0.03649 
## F-statistic: 3.575 on 3 and 201 DF,  p-value: 0.01493

Adjusted R2-value of the model was 0.03649, and multiple R2 was 0.05066. The negative adjusted R2-value indicates that the amount of variation accounted for by Cover crops is negligible. In practice, cover crops explain no variation in the data.

The effect sizes (marginal means) of % NH3-loss for each cover category can be seen below.

9Model 7-8 - Air temperature and rainfall

AirTemperature over the duration of the experiment (Celcius) and mean rainfall over the experimental period (mm) were analysed together, to explore interactions between the two variables. The frequency distributions of the two variables can be seen below.

There are 154 observations in the analysis.

# Model
mTR_lm <- lm(sqrt(NH3loss+4) ~ AirTemperature * Rain.mean + offset(M1_fit), 
             weights=Weight, 
             data=NH3predict_m7)

The model assumptions (heteroscedasticity and normality of errors):

There was a marginally non-significant effect of interaction (F(1, 150) =6, p = 0.015), and a marginally non-significant main effect of mean rainfall (F(1, 150) =0.55, p = 0.457). There was no main effect of air temperature (F(1, 150) =5.7, p = 0.018).

Full model output can be seen below:

## 
## Call:
## lm(formula = sqrt(NH3loss + 4) ~ AirTemperature * Rain.mean + 
##     offset(M1_fit), data = NH3predict_m7, weights = Weight)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3688 -0.8447 -0.1843  0.7469  4.7047 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)  
## (Intercept)              -0.066938   0.211120  -0.317   0.7516  
## AirTemperature            0.009843   0.014395   0.684   0.4952  
## Rain.mean                -0.179586   0.072468  -2.478   0.0143 *
## AirTemperature:Rain.mean  0.009858   0.004026   2.449   0.0155 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.263 on 150 degrees of freedom
## Multiple R-squared:  0.07554,    Adjusted R-squared:  0.05705 
## F-statistic: 4.085 on 3 and 150 DF,  p-value: 0.008026

Air temperature and rainfall, and the interaction between the two explain 5.7 % of the variation in the data.

## Using data NH3predict_m7 from global environment. This could cause
## incorrect results if NH3predict_m7 has been altered since the model was
## fit. You can manually provide the data to the "data =" argument.
## Outcome is based on a total of 1 exposures.

## Using data NH3predict_m7 from global environment. This could cause
## incorrect results if NH3predict_m7 has been altered since the model was
## fit. You can manually provide the data to the "data =" argument.
## Outcome is based on a total of 1 exposures.

## Using data NH3predict_m7 from global environment. This could cause
## incorrect results if NH3predict_m7 has been altered since the model was
## fit. You can manually provide the data to the "data =" argument.
## Loading required namespace: broom.mixed

## JOHNSON-NEYMAN INTERVAL 
## 
## When Rain.mean is OUTSIDE the interval [-14.37, 1.46], the slope of
## AirTemperature is p < .05.
## 
## Note: The range of observed values of Rain.mean is [0.00, 12.00]

## Using data NH3predict_m7 from global environment. This could cause
## incorrect results if NH3predict_m7 has been altered since the model was
## fit. You can manually provide the data to the "data =" argument.

## JOHNSON-NEYMAN INTERVAL 
## 
## When AirTemperature is OUTSIDE the interval [9.12, 37.96], the slope of
## Rain.mean is p < .05.
## 
## Note: The range of observed values of AirTemperature is [-5.00, 27.70]

Air temperature as a separate model

Considering air temperature separately, there is a statistically significant effect of air temperature (F(1, 159) =6.47, p = 0.012). Adjusted R2-value of the model was 0.03307, and multiple R2 was 0.03912. This means air temperature explains 3.307 % of the variation.

There are 161 observations in the analysis.

The full model output with the slopes and their standard error can be seen below:

## 
## Call:
## lm(formula = sqrt(NH3loss + 4) ~ AirTemperature + offset(M1_fit), 
##     data = NH3predict_m7_2, weights = Weight)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3336 -0.9530 -0.1121  0.7980  4.6170 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    -0.40256    0.16559  -2.431   0.0162 *
## AirTemperature  0.03044    0.01196   2.544   0.0119 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.26 on 159 degrees of freedom
## Multiple R-squared:  0.03912,    Adjusted R-squared:  0.03307 
## F-statistic: 6.473 on 1 and 159 DF,  p-value: 0.01191

Air temperature over first 7 days

There is no statistically significant effect of air temperature over the first 7 days of the experiment (F(1, 18) =1.77, p = 0.2). Adjusted R2-value of the model was 0.039, and multiple R2 was 0.08958. This means air temperature over the first 7 days explains 3.9 % of the variation.

There are 20 observations in the analysis.

The full model output with the slopes and their standard error can be seen below:

## 
## Call:
## lm(formula = sqrt(NH3loss + 4) ~ Temp.7d + offset(M1_fit), data = NH3predict_mA7, 
##     weights = Weight)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.5171 -1.1222 -0.3516  1.0465  4.5712 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.48292    0.53642   0.900     0.38
## Temp.7d     -0.09001    0.06764  -1.331     0.20
## 
## Residual standard error: 1.749 on 18 degrees of freedom
## Multiple R-squared:  0.08958,    Adjusted R-squared:  0.039 
## F-statistic: 1.771 on 1 and 18 DF,  p-value: 0.1999

Mean rainfall as a separate model

Considering mean rainfall separately, there is no statistically significant effect of rain (F(1, 185) =0.29, p = 0.591). Adjusted R2-value of the model was -0.00383, and multiple R2 was 0.00157. Thus, rainfall doesn’t meaningfully explain any variation in the data.

There are 187 observations in the analysis.

The full model output with the slopes and their standard error can be seen below:

## 
## Call:
## lm(formula = sqrt(NH3loss + 4) ~ Rain.mean + offset(M1_fit), 
##     data = NH3predict_m8, weights = Weight)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4817 -0.8414 -0.1172  0.7477  5.0631 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.01874    0.08881  -0.211    0.833
## Rain.mean   -0.01963    0.03643  -0.539    0.591
## 
## Residual standard error: 1.259 on 185 degrees of freedom
## Multiple R-squared:  0.001567,   Adjusted R-squared:  -0.00383 
## F-statistic: 0.2904 on 1 and 185 DF,  p-value: 0.5906

Rainfall over first 72 h

There is no statistically significant effect of rainfall during the first 72h of the experiment (F(1, 127) =0.55, p = 0.46). Adjusted R2-value of the model was -0.00353, and multiple R2 was 0.00431. Thus, rainfall during the first 72h doesn’t meaningfully explain any variation in the data.

There are 129 observations in the analysis.

The full model output with the slopes and their standard error can be seen below:

## 
## Call:
## lm(formula = sqrt(NH3loss + 4) ~ Rain.72h + offset(M1_fit), data = NH3predict_R8, 
##     weights = Weight)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3564 -0.9793 -0.1928  0.7826  5.0386 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.01284    0.09684  -0.133    0.895
## Rain.72h    -0.01475    0.01991  -0.741    0.460
## 
## Residual standard error: 1.328 on 127 degrees of freedom
## Multiple R-squared:  0.004306,   Adjusted R-squared:  -0.003535 
## F-statistic: 0.5492 on 1 and 127 DF,  p-value: 0.46

Figures below show a linear approximation of the effect of air temperature (Celcius) and mean rainfall (mm) on NH3-emissions (linear slope back-transformed from the square-root scale).

Rainfall as a non-linear effect

We fitted an exponential function with parameters a and b to rainfall data. The equation is as follows:

\(f(Rain.mean) = a * e^{(b*Rain.mean)}\)

The response variable is presented in square-root(+4) scale to make the data comparable for further analysis.

To test the fit and to estimate parameters a and b, we fitted a non-linear least-squares model.

## 
## Formula: NH3loss ~ a * exp(b * (Rain.mean))
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a 13.20159    1.26103  10.469   <2e-16 ***
## b  0.01001    0.03539   0.283    0.778    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.88 on 152 degrees of freedom
## 
## Number of iterations to convergence: 14 
## Achieved convergence tolerance: 1.49e-08

The fitted relationship is as below:

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

The RMSE (root mean square error) of the model was 12.7975579, which is on the same scale as mean NH3loss values in the data (13.4038961) indicating a poor fit of the model. (Error equal to the estimated value.)

We calculated new rain estimates on the exponential scale, with the estimated parameters taken into account, and ran the analysis with air temperature again, using the exponential Rain.mean values calculated using the coefficients (parameter estimates) from the non-linear least-squares fit above.

# Transforming the Rain values to exp-scale
NH3predict_m7$Rain.exp <- coef(MR_nls)[[1]]*exp(coef(MR_nls)[[2]]*NH3predict_m7$Rain.mean)

## 
## Call:
## lm(formula = sqrt(NH3loss + 4) ~ AirTemperature * Rain.exp + 
##     offset(M1_fit), data = NH3predict_m7, weights = Weight)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3699 -0.8446 -0.1871  0.7469  4.7070 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)  
## (Intercept)             17.05144    7.08396   2.407   0.0173 *
## AirTemperature          -0.92197    0.38988  -2.365   0.0193 *
## Rain.exp                -1.29721    0.52690  -2.462   0.0150 *
## AirTemperature:Rain.exp  0.07062    0.02894   2.440   0.0159 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.263 on 150 degrees of freedom
## Multiple R-squared:  0.0751, Adjusted R-squared:  0.0566 
## F-statistic:  4.06 on 3 and 150 DF,  p-value: 0.008292

When the interaction is removed, we end in a model showing very similar results as the linear model above, with no statistically significant effect of rainfall (F(1, 151) =0.51, p = 0.47), nor air temperature (F(1, 151) =5.52, p = 0.018), but a statistically significant slope for air temperature (see below).

## 
## Call:
## lm(formula = sqrt(NH3loss + 4) ~ AirTemperature + Rain.exp + 
##     offset(M1_fit), data = NH3predict_m7, weights = Weight)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3589 -0.9590 -0.1725  0.7856  4.6484 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)  
## (Intercept)     2.26437    3.72808   0.607   0.5445  
## AirTemperature  0.02889    0.01233   2.343   0.0204 *
## Rain.exp       -0.19804    0.27777  -0.713   0.4770  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.284 on 151 degrees of freedom
## Multiple R-squared:  0.03839,    Adjusted R-squared:  0.02565 
## F-statistic: 3.014 on 2 and 151 DF,  p-value: 0.05205

When non-linear rain is analysed separately, there is no effect of rain on NH3loss (F(1, 185) =0.26, p = 0.609). See the slope estimate from below:

## 
## Call:
## lm(formula = sqrt(NH3loss + 4) ~ Rain.exp + offset(M1_fit), data = NH3predict_m8, 
##     weights = Weight)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4816 -0.8415 -0.1242  0.7494  5.0625 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   1.7636     3.5368   0.499    0.619
## Rain.exp     -0.1352     0.2638  -0.512    0.609
## 
## Residual standard error: 1.259 on 185 degrees of freedom
## Multiple R-squared:  0.001417,   Adjusted R-squared:  -0.003981 
## F-statistic: 0.2624 on 1 and 185 DF,  p-value: 0.6091

Model 9 -Irrigation

Irrigation was modeled as a binomial variable due to heavy zero-inflation in the data.

The number of observations in total in each category can be seen below:

# Model 
m9_lm <- lm(sqrt(NH3loss+4) ~ irrigation_cat + offset(M1_fit), 
            weights=Weight,
            data=NH3predict_m9)

The check for normality and heteroscedasticity of errors can be seen below:

There was no effect of irrigation on NH3-loss (F(1, 162) =0.02, p = 0.886).

## 
## Call:
## lm(formula = sqrt(NH3loss + 4) ~ irrigation_cat + offset(M1_fit), 
##     data = NH3predict_m9, weights = Weight)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3338 -0.9997 -0.1092  0.8246  5.0354 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)
## (Intercept)     -0.03160    0.07984  -0.396    0.693
## irrigation_cat1 -0.05269    0.36787  -0.143    0.886
## 
## Residual standard error: 1.331 on 162 degrees of freedom
## Multiple R-squared:  0.0001266,  Adjusted R-squared:  -0.006045 
## F-statistic: 0.02052 on 1 and 162 DF,  p-value: 0.8863

Adjusted R2-value of the model was -0.00605, and multiple R2 was 1.3^{-4}. The negative adjusted R2-value indicates that the amount of variation accounted for by irrigation is negligible. In practice, it explain no variation in the data.

The marginal means for NH3 loss in no-irrigation and irrigation categories can be seen below.

Model 10 - Water input

Daily water input (mm) consists of irrigation and rainfall.

#  Model
m10_lm <- lm(sqrt(NH3loss+4) ~ WaterPerDay + offset(M1_fit), 
             weights=Weight,
             data=NH3predict_m10)

Heteroscedasticity and normality of erros were as below:

There was no effect of water input on NH3-loss (F(1, 189) =0.41, p = 0.525).

## 
## Call:
## lm(formula = sqrt(NH3loss + 4) ~ WaterPerDay + offset(M1_fit), 
##     data = NH3predict_m10, weights = Weight)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4639 -0.8855 -0.0893  0.7717  5.0934 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.03447    0.08695  -0.396    0.692
## WaterPerDay -0.02119    0.03326  -0.637    0.525
## 
## Residual standard error: 1.295 on 189 degrees of freedom
## Multiple R-squared:  0.002143,   Adjusted R-squared:  -0.003136 
## F-statistic: 0.4059 on 1 and 189 DF,  p-value: 0.5248

Adjusted R2-value of the model was -0.00314, and multiple R2 was 0.00214. The negative adjusted R2-value indicates that the amount of variation accounted for by daily water input is negligible. In practice, it explain no variation in the data.

Figure below shows a linear approximation of the effect of daily water input (mm) on NH3-emissions (linear slope back-transformed from the square-root scale).

Model 11 - Experiment duration

Duration of the experiment (days).

#  Model
m11_lm <- lm(sqrt(NH3loss+4) ~ Duration + offset(M1_fit), 
             weights=Weight,
             data=NH3predict_m11)

Heteroscedasticity and normality of erros were as below:

There was no effect of duration of the experiment on NH3-loss (F(1, 193) =0.12, p = 0.727).

## 
## Call:
## lm(formula = sqrt(NH3loss + 4) ~ Duration + offset(M1_fit), data = NH3predict_m11, 
##     weights = Weight)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4700 -0.8108 -0.0885  0.7788  5.0518 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.070463   0.112494  -0.626    0.532
## Duration     0.001077   0.003079   0.350    0.727
## 
## Residual standard error: 1.242 on 193 degrees of freedom
## Multiple R-squared:  0.0006338,  Adjusted R-squared:  -0.004544 
## F-statistic: 0.1224 on 1 and 193 DF,  p-value: 0.7268

Adjusted R2-value of the model was -0.00454, and multiple R2 was 6.3^{-4}. The negative adjusted R2-value indicates that the amount of variation accounted for by duration of the experiment is negligible. In practice, it explain no variation in the data.

Figure below shows a linear approximation of the effect of daily water input (mm) on NH3-emissions (linear slope back-transformed from the square-root scale).

Combined model with all variables of effect except application rate

Model 12 - Base model, pH, and temperature

Total number of observations 155

mALL_lmer <- lmer(sqrt(NH3loss+4) ~ Fertiliser + Application + 
                    pH_cat + Fertiliser:pH_cat + AirTemperature +
                    (1|Ref),
                   weights=Weight,
                   control=lmerControl(optimizer="bobyqa"),
                   data=NH3predict_mALL)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## sqrt(NH3loss + 4) ~ Fertiliser + Application + pH_cat + Fertiliser:pH_cat +  
##     AirTemperature + (1 | Ref)
##    Data: NH3predict_mALL
## Weights: Weight
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 468.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1370 -0.5491 -0.0660  0.5269  3.2042 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Ref      (Intercept) 0.8238   0.9077  
##  Residual             1.6150   1.2708  
## Number of obs: 155, groups:  Ref, 27
## 
## Fixed effects:
##                                    Estimate Std. Error        df t value
## (Intercept)                         3.85697    0.33669  55.86542  11.456
## FertiliserUAN                      -0.71093    0.52992 143.78255  -1.342
## Fertiliserammonium+1               -1.45739    0.63015 134.14941  -2.313
## Fertiliserammonium+2               -0.82239    0.77588 128.01895  -1.060
## ApplicationIncorporated            -1.34907    0.66476 141.75588  -2.029
## ApplicationLiquid                  -0.50287    0.37641 140.28090  -1.336
## pH_catchalky                       -0.65003    0.27315 141.11050  -2.380
## AirTemperature                      0.05391    0.01723 143.00870   3.128
## FertiliserUAN:pH_catchalky          1.03221    0.53827 141.18035   1.918
## Fertiliserammonium+1:pH_catchalky   1.25151    0.66947 141.12157   1.869
## Fertiliserammonium+2:pH_catchalky   0.66515    0.80195 127.77651   0.829
##                                   Pr(>|t|)    
## (Intercept)                       2.74e-16 ***
## FertiliserUAN                      0.18185    
## Fertiliserammonium+1               0.02226 *  
## Fertiliserammonium+2               0.29116    
## ApplicationIncorporated            0.04429 *  
## ApplicationLiquid                  0.18373    
## pH_catchalky                       0.01866 *  
## AirTemperature                     0.00213 ** 
## FertiliserUAN:pH_catchalky         0.05718 .  
## Fertiliserammonium+1:pH_catchalky  0.06364 .  
## Fertiliserammonium+2:pH_catchalky  0.40842    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) FrtUAN Frtl+1 Frtl+2 ApplcI ApplcL pH_ctc ArTmpr FUAN:H
## FertilsrUAN -0.133                                                        
## Frtlsrmmn+1 -0.253  0.035                                                 
## Frtlsrmmn+2 -0.086  0.239  0.077                                          
## ApplctnIncr  0.034 -0.001 -0.010 -0.016                                   
## ApplictnLqd -0.182 -0.459  0.141 -0.324 -0.003                            
## pH_catchlky -0.258  0.182  0.262  0.264 -0.035  0.040                     
## AirTempertr -0.639 -0.001  0.079  0.012 -0.107  0.102 -0.136              
## FrtlsUAN:H_  0.214 -0.728 -0.071 -0.021  0.009 -0.059 -0.289 -0.034       
## Frtlsr+1:H_  0.247 -0.085 -0.897 -0.083  0.022  0.011 -0.303 -0.116  0.133
## Frtlsr+2:H_  0.125 -0.062 -0.094 -0.823  0.018 -0.020 -0.333 -0.012  0.100
##             F+1:H_
## FertilsrUAN       
## Frtlsrmmn+1       
## Frtlsrmmn+2       
## ApplctnIncr       
## ApplictnLqd       
## pH_catchlky       
## AirTempertr       
## FrtlsUAN:H_       
## Frtlsr+1:H_       
## Frtlsr+2:H_  0.110
##             R2m       R2c
## [1,] 0.09943048 0.4036482

Model 11 explains 40.36 % of the variation in the data (conditional R2), with the fixed effects explaining 9.94 % (marginal R2).

The same effect sizes (and their 95% confidence interval) are shown below:

Combined model with all variables of effect, including application rate

Model 13 - Base model, pH, rainfall and temperature

#
# Remove NA's from the dataset
NH3predict_mALL2 <- NH3predict[which(is.na(NH3predict$pH_cat)==FALSE &  is.na(NH3predict$AirTemperature)==FALSE &  is.na(NH3predict$Rainfall)==FALSE),]

# Transforming the Rain values to exp-scale
NH3predict_mALL2$Rain.exp <- coef(MR_nls)[[1]]*exp(coef(MR_nls)[[2]]*NH3predict_mALL2$Rain.mean)


mALL_lmer2 <- lmer(sqrt(NH3loss+4) ~ Fertiliser + Application + 
                    pH_cat + Fertiliser:pH_cat + AirTemperature + Rain.exp+
                     AirTemperature:Rain.exp+
                    (1|Ref),
                   weights=Weight,
                   control=lmerControl(optimizer="bobyqa"),
                   data=NH3predict_mALL2)

Total observations of150

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## sqrt(NH3loss + 4) ~ Fertiliser + Application + pH_cat + Fertiliser:pH_cat +  
##     AirTemperature + Rain.exp + AirTemperature:Rain.exp + (1 |      Ref)
##    Data: NH3predict_mALL2
## Weights: Weight
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 449.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1761 -0.5156 -0.0839  0.5212  3.3611 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Ref      (Intercept) 0.6391   0.7994  
##  Residual             1.6248   1.2747  
## Number of obs: 150, groups:  Ref, 26
## 
## Fixed effects:
##                                    Estimate Std. Error        df t value
## (Intercept)                        33.53193   10.76764 107.60554   3.114
## FertiliserUAN                      -0.66367    0.52532 136.87630  -1.263
## Fertiliserammonium+1               -0.96185    0.65110 133.07711  -1.477
## Fertiliserammonium+2               -0.92596    0.77743 122.11360  -1.191
## ApplicationIncorporated            -1.69656    0.68125 135.98575  -2.490
## ApplicationLiquid                  -0.47626    0.37457 135.49051  -1.271
## pH_catchalky                       -0.62740    0.27149 130.40540  -2.311
## AirTemperature                     -1.48557    0.51222 131.41267  -2.900
## Rain.exp                           -2.21547    0.80319 106.74025  -2.758
## FertiliserUAN:pH_catchalky          0.98737    0.53079 129.37002   1.860
## Fertiliserammonium+1:pH_catchalky   0.80759    0.68775 136.50481   1.174
## Fertiliserammonium+2:pH_catchalky   0.71438    0.80357 121.80973   0.889
## AirTemperature:Rain.exp             0.11490    0.03823 130.92076   3.006
##                                   Pr(>|t|)   
## (Intercept)                        0.00236 **
## FertiliserUAN                      0.20861   
## Fertiliserammonium+1               0.14197   
## Fertiliserammonium+2               0.23594   
## ApplicationIncorporated            0.01397 * 
## ApplicationLiquid                  0.20574   
## pH_catchalky                       0.02240 * 
## AirTemperature                     0.00437 **
## Rain.exp                           0.00684 **
## FertiliserUAN:pH_catchalky         0.06513 . 
## Fertiliserammonium+1:pH_catchalky  0.24234   
## Fertiliserammonium+2:pH_catchalky  0.37575   
## AirTemperature:Rain.exp            0.00318 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

Model 12 explains 37.4 % of the variation in the data (conditional R2), with the fixed effects explaining 12.77 % (marginal R2).

The same effect sizes (and their 95% confidence interval) are shown below:

References

Astivia, Oscar L. Olvera and Zumbo, Bruno D. (2019) “Heteroskedasticity in Multiple Regression Analysis: What it is, How to Detect it and How to Solve it with Applications in R and SPSS,” Practical Assessment, Research, and Evaluation: Vol. 24 , Article 1. DOI: https://doi.org/10.7275/q5xr-fr95

Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1-48. https://doi:10.18637/jss.v067.i01

Nakagawa, S, Schielzeth, H. (2013). A general and simple method for obtaining R^2 from Generalized Linear Mixed-effects Models. Methods in Ecology and Evolution 4: 133-142. https://doi.org/10.1111/j.2041-210x.2012.00261.x